#include "KDM1PrimaryGeneratorAction.hh"
#include "KDM1Globals.hh"
#include "G4ParticleTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4ThreeVector.hh"
#include "G4ios.hh"
#include "G4HEPEvtInterface.hh"
#include "G4Event.hh"
#include "G4PrimaryVertex.hh"

KDM1PrimaryGeneratorAction::KDM1PrimaryGeneratorAction()
{
	G4Event* tempEvent = new G4Event();
	G4HEPEvtInterface* eventInterface = new G4HEPEvtInterface(eventInterfaceFilePath);
	numOfVertices = kaonCount*3;
	vertexList = new G4PrimaryVertex*[numOfVertices];
	energyList = new double[numOfVertices];
	std::ifstream file;
	file.open(kaonPositionFilePath);
	double c = 299792458*m/s;
	double kLongMass = 497.648*MeV/c/c;
	double e = 2.71828183;
	double tao = 51.16*ns;
	for(int i = 0;i < kaonCount;i++)
	{
		float x = 0;
		float y = 0;
		float z = 0;
		float vx = 0;
		float vy = 0;
		float vz = 0;
		int counter = 0;
		while(counter < 7)
		{
			char* number = new char[64];
			file.getline(number,64,' ');
			while(!(number[0])) file.getline(number,64,' ');
			float numberf = atof(number);
			switch(counter)
			{
				case 0: x = numberf*mm;
					break;
				case 1: y = numberf*mm;
					break;
				case 2: z = numberf*mm-2.75*m;
					break;
				case 3: vx = (numberf*GeV/c)/sqrt(kLongMass*kLongMass+pow(numberf*GeV/c/c,2));
					break;
				case 4: vy = (numberf*GeV/c)/sqrt(kLongMass*kLongMass+pow(numberf*GeV/c/c,2));
					break;
				case 5: vz = (numberf*GeV/c)/sqrt(kLongMass*kLongMass+pow(numberf*GeV/c/c,2));
					break;
				case 6: energyList[i] = numberf*GeV;
					break;
			}
			counter++;
			delete number;
		}
		double maxPropagationTime = 1.6*m/vz;
		double propagationTime = tao*log(1/(1-(1-pow(e,-maxPropagationTime/tao))*(i+.5)/kaonCount));
		x += vx*propagationTime;
		y += vy*propagationTime;
		z += vz*propagationTime;
		eventInterface->GeneratePrimaryVertex(tempEvent);
		G4PrimaryVertex* tempVertex = tempEvent->GetPrimaryVertex(i);
		for(int j = 0;j < 3;j++)
		{
			vertexList[3*i+j] = new G4PrimaryVertex(G4ThreeVector(x,y,z),propagationTime);
			G4PrimaryParticle* tempParticle = tempVertex->GetPrimary(j);
			G4PrimaryParticle* tempParticleNew = new G4PrimaryParticle(tempParticle->GetG4code());
			double px = tempParticle->GetPx();
			double py = tempParticle->GetPy();
			double pz = tempParticle->GetPz();
			tempParticleNew->SetMomentum(px,py,pz);
			vertexList[3*i+j]->SetPrimary(tempParticleNew);
		}
	}

	delete eventInterface;
	delete tempEvent;
}

KDM1PrimaryGeneratorAction::~KDM1PrimaryGeneratorAction()
{
	delete vertexList;
}

void KDM1PrimaryGeneratorAction::GeneratePrimaries(G4Event* event)
{
	eventNum++;
	event->AddPrimaryVertex(vertexList[eventNum]);
}
